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ABSTRACT 

We have recently interpreted the source MAGIC J0616+225 as a result of de- 
layed TeV emission of cosmic -rays diffusing from IC 443 and interacting with a 
cloud in the foreground of the remnant. This model was used to make predictions 
for future observations, especially those to be made with the Fermi satellite. Just re- 
cently, AGILE, Fermi, and VERITAS have released new results of their observations 
of IC 443. In this work, we compare them with the predictions of our model, explor- 
ing the GeV to TeV connection in this region of space. We use Fermi data to consider 
the possibility of constraining the cosmic -ray diffusion features of the environment. 
We analyze the cosmic -ray distributions, their interactions, and a possible detection 
of the SNR environment in the neutrino channel. 

Key words: SNR (individual IC 443), y-rays: observations, y-rays: theory 



1 INTRODUCTION 

It is commonly accepted that supernova remnants (SNR) are one of the most probable scenarios 
of leptonic and hadronic cosmic-ray (CR) acceleration. The particle acceleration mechanism in 
individual SNRs is usually assumed to be diffusive shock acceleration, which naturally leads to 
a power-law population of relativistic particles. In the standard version of this mechanism (e.g. 
Bell 1978), particles are scattered by magnetohydrodynamic waves repeatedly through the shock 
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front. Electrons suffer synchrotron losses, producing the non-thermal emission from radio to X- 
rays usually seen in shell-type SNRs. The maximum energy achieved depends on the shock speed 
and age as well as on any competing loss processes. In young SNRs, electrons can easily reach 
energies in excess of 1 TeV, and they produce X-rays. Non-thermal X-ray emission associated with 
shock acceleration has been clearly observed in many SNRs. But in order to have an observational 
confirmation of protons and other nuclei being accelerated, particularly, in order to be able to 
distinguish this from leptonic emission, one should try and isolate the multi-messenger effects 
of the secondary particles produced when the accelerated hadrons interact in nearby molecular 
clouds through pp collisions. These ideas go back, for instance, to the works by Dogel & Sharov 
1990, Naito & Takahara 1994; Drury 1994; Sturner et al. 1997; Gaisser et al. 1998; Baring et al. 
1999, among others. In fact, as early as 1979, Montmerle suggested that SNRs within OB stellar 
associations, i.e. star forming regions with plenty of molecular gas, could generate observable y- 
ray sources. A molecular cloud being illuminated by particles that escaped from a nearby SNR 
could then act as a target for pp interactions, greatly enhancing the y-ray emission (see, e.g., the 
recent works by Gabici et al. 2007, 2009; Casanova et al. 2009; Rodriguez-Marrero et al. 2008). 
As an spinoff, observing y-rays from clouds nearby SNRs, can feedback on our knowledge of 
the diffusion characteristics of the environment. As has been emphasized by Aharonian & Atoyan 
(1996), the observed y-ray s can have a significantly different spectrum from that expected from 
the primary particle population at the immediate vicinity of source (the SNR shock). For instance, 
a standard diffusion coefficients 8 ~ 0.3 - 0.6 can explain y-ray spectra as steep as Y ~ 2.3 - 2.6 in 
sources with particles accelerated to a power-law J P (E P ) oc E~ 2 if the target that is illuminated by 
the 7r°-decays is at sufficient distance from the accelerator. Measuring y-ray emission around SNRs 
would then allow to acquire knowledge of the diffusion environment in which the CRs propagate, 
at several kpc from Earth. 

Of all SNRs that were found to be positionally coincident with y-ray sources in the MeV 
range in the EGRET era, IC 443 was one of the most appealing for subsequent observations with 
higher sensitivity instruments (see the case-by-case study by Torres et al. 2003). It was, perhaps 
with W28, the only case in which the molecular environment -as mapped for instance with CO 
observations- showed a peak in density close by, but separated in sky projection, from the SNR 
center. This would allow distinguishing, in case the y-ray emission observed would be hadronically 
produced, possible cosmic-ray diffusion effects. Along the last year, several new observations of 
the IC 443 environment have been made, and in this work, we consider these in the setting of a 
theoretical model in which CRs from the SNR IC 443 are diffusing away from it and interacting 
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with clouds nearby. This model was originally put forward by Aharonian & Atoyan (1996), and 
Torres et al. (2008), referred to as Paper I in this work, studied this model for IC 443 prior to the 
new wealth of data we can now consider. 

2 HIGH AND VERY HIGH-ENERGY OBSERVATIONS 

2.1 Earlier EGRET and MAGIC observations 

MAGIC observations towards IC 443 yielded the detection of J06 16+225 nearby, but displaced 
from the center of the SNR IC 443, with centroid located at (RA,DEC) 72 ooo=(06 h 16 m 43 s , +22°31' 
48"), +0.025° to , + 0.017° ys (Albert et al. 2007). No extension nor any variability was claimed 
in the y-ray data. Albert et al. (2007) showed that the MAGIC source is located at the position 
of a giant cloud in front of the SNR. A simple power law was fitted to the measured spectral 
points: dN 7 /(dAdtdE) = (1.0 + 0.2 stat + 0.35,„) x 10" 11 (£/0.4TeV)" 31±0 - 3 "" f±0 - 2 '" cnrV'TeV -1 . 
The integral flux of MAGIC J0616+225 above 100 GeV is about 6.5% of the Crab Nebula. The 
EGRET flux of the source 3EG J06 17+2238, which is positionally correlated with the SNR IC 443, 
is (51.4+3.5) xlO -8 ph cm~ 2 s^ 1 , and it presents a photon spectral index of 2.01+0.06 (Hartman et 
al. 1999). The EGRET source was classified as non-variable by Torres et al. (2001) and Nolan et 
al. (2003). An independent analysis of GeV photons measured by EGRET resulted in the source 
GeV J0617+2237 (Lamb & Macomb 1997), also at the same location of 3EG J0617+2238, the 
centroid of which is at the center of the SNR shell. 

2.2 Recent TeV observations 

Recently, the Very Energetic Radiation Imaging Telescope Array System (VERITAS) presented 
further observations towards IC 443 (Acciari et al. 2009). Regarding the position of the centroid, 
it was found to be at (RA,DEC) 72 ooo=(06 h 16 m 51 s ,+22°30' 11"), ±0.03° te , + 0.08° }>? thus, consistent 
with that of MAGIC. Evidence that the very-high-energy (VHE, E > 100 GeV) y-ray emission is 
extended was also found. The extension derived was 0.16° + 0.03° te; + 0.04° . The VHE spectrum 
is well fit by a power law (dN/dE = N x (£/TeVr r ) with a photon index of 2.99 + 0.38^ + 0.3 sy . s 
and an integral flux above 300 GeV of (4.63 + 0.90 stat + 0.93, y ,) x 10~ 12 cnr 2 s _1 . Thus, as we will 
graphically see below, the spectral determination is consistent with the MAGIC measurements, 
both present a steep slope, with VERITAS finding a slight overall increase in the flux level. No 
variability of the y-ray emission was claimed by VERITAS either. 
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2.3 Recent GeV observations 

AGILE results on IC 443 has been recently reported too (Tavani et al. 2010). AGILE discovered 
a distinct pattern of diffuse emission in the energy range 100 MeV-3 GeV coming from the SNR, 
with a prominent maximum localized in the Northeastern shell, dislocated (as it was the case with 
EGRET) with the MAGIC/VERITAS sources. The latter is -0.4° apart from the maximum of the 
AGILE emission (which in turn is also away from the nearby PWN, discussed below). Finally, 
Fermi has also recently presented an analysis of its first 1 1 months of observations towards the 
region of SNR IC 443 (Abdo et el. 2010). These results enhance, given the better instrument 
sensitivity, those obtained by AGILE. Thus, we focus on Fermi measurements when analyzing 
GeV results. The source was detected in a broad range of energies, from 200 MeV up to 50 GeV, 
with a SED that rolls over at about 3 GeV to seemingly match in slope the one that is found at the 
highest energies; i.e., it can be represented, for instance, with a broken power law with slopes of 
1.93 + 0.03 and 2.56 + 0.1 1 and with a break at 3.25 + 0.6 GeV. This is one important difference 
with EGRET data, which SED did not allow to suspect neither that the emission would maintain 
a hard spectrum up to such tens-of-GeV energies nor the existence of a roll over in the spectrum 
at the energies found. The flux above 200 MeV resulted to be (28.5+0.7) xlO -8 ph cnT 2 s _1 what 
allowed for a very significant detection in Fermi. The centroid of the emission is consistent with 
that of EGRET 3EG J06 17+2238. 

2.4 Relative localization of sources 

Abdo et al. (2010) report that the centroid of the Fermi emission is displaced more than 5 x 
^ r 8 ror (MAGIC error) from that of MAGIC (J0610+225), and more than 1.5 x 0%™ (VERITAS 
error) from that of the VERITAS source. These numbers are obtained assuming that the systematic 
and statistical errors in localization add up in quadrature, and considering the worse error of each 
of the pairs of measurements (Fermi-MAGIC, Fermi- VERITAS), which in both cases correspond 
to the IACTs. The significance of the separation greatly improves when a) the best measured 
position is considered (i.e., the error by Fermi), for which both pairs of measurements are about 
5cr away, and/or b) when statistical errors only are considered for VERITAS (the systematic errors 
in this latter measurement is about a factor of 3 larger than the statistics and significantly different 
from all others, but of course, one can not necessarily assume it to approach the detection in the 
direction of the Fermi source). Thus, albeit current measurements are not conclusive about energy 
dependent morphology, they are consistent with it: Abdo et al. (2010) report that the centroid of 
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the emission moves (but still not significantly in Fermi data: only at ~ 1.5<x) towards that of the 
VERITAS source as the energy band changes from 1-5 GeV to 5-50 GeV It might be that the 
angular resolution and/or sensitvity and/or the separation of the real molecular mass distribution 
on sky projection are not enough to distinguish the difference when such nearby energy ranges are 
considered. New measurements from MAGIC (using the just-obtained stereoscopic capability of 
the array) could provide continuous coverage from 50 GeV up. 

2.5 APWN? 

In all energy bands, the centroid of the correspondingly detected sources is inconsistent with the 
pulsar wind nebula (and the putative pulsar) CXOU J06 1705. 3 +222 127, discovered by Olbert et 
al. (2001), and lying nearby. Both the 3EG and the GeV source in the catalogs of Hartman et al. 
(1999) and Lamb & Macomb (1997), which are co-spatial, are inconsistent with the PWN location. 
Similarly, using the position of the PWN and the Fermi source one sees that they are separated by 
0.26°, or about llcr away from the localization of the Fermi peak. Also at higher energies, the 
y-ray emission observed by VERITAS and MAGIC is offset from the location of the PWN by 
10-20 arcmin. This latter fact could be understood in case the PWN is a y-ray emitter, it would be 
similar to the case of HESS J1825-137 (Aharonian et al. 2006) or HESS J1908+063 (Aharonian 
et al. 2009), where similar offsets were found, see also Abdo et al. (2010b). The emission could 
be consistent with a scenario in which the VHE emission arises from inverse Compton scattering 
off electrons accelerated early in the PWN's life. However, if one would assume that the PWN 
CXOU J06 1705. 3 +222 127 is producing the emission (note that pulsed radiation from this object 
has not been found at any frequency), the highest energy TeV-band radiation should peak there (it 
could be extended, but due to losses, the higher the energy, the more peaked towards the PWN the 
emission will be) and the GeV radiation should then be unresolved pulsar emission, it should also 
peak there and be pulsed (see Bartko & Bednarek 2008). Then, based on FermifVERlTAS data we 
can safely entertain that the GeV and TeV emissions detected do not originate in the PWN, what 
we explore here further. 

3 THE COSMIC-RAY DIFFUSION MODEL 

As a first approach to the modeling of the GeV to TeV SED, Abdo et al. (2010) assumed that a 
single proton spectrum was directly interacting with the whole molecular mass found in the IC 443 
environment (see Torres et al. 2003). Given that the location of sources at different energies change, 
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and that target material for accelerated cosmic-rays are also found to be at different positions, 
this approach is just a crude approximation to the need of considering cosmic-ray diffusion. This 
kind of model was introduced in Paper I, and we refer the reader there for details: basically, we 
computed the spectrum of y-rays generated through 7r°-decay at a source of proton density n p (see 
e.g., Torres 2004 or Domingo-Santamaria & Torres 2005 for the formulae we used); solving for 
the cosmic ray spectrum at each distance of the SNR, neglecting temporal or spatial effects (non- 
uniform densities) within the molecular cloud itself. The spectrum of y-rays generated through 
7r°-decay at a source of proton density n p is 

F y(E 7 ) = 2 [ {F n {En)l J El -ml) dE n , (1) 
where the minimum pion energy is E™ m {E y ) = E 7 + ml/4E y , and 

F n {E n ) = 4nn p \ " J p (E)(do- n (E n , E p )/dE n ) dE p . (2) 

Here, dcr„(E K , E p )/dE K is the differential cross-section for the production of 7T°-mesons of energy 
E n by a proton of energy E p in a pp collision. We analyze below the influence upon the results of 
different parameterizations of this cross section. We include cosmic-rays and target nuclei heavier 
than the proton throughout this paper. We adopt here the approximation in which these can be 
accounted for by multiplying a nuclear factor (adopted as 1.5) to the total flux, without changing 
the cosmic -ray proton spectrum (Gaisser & Schaefer 1992). 

We have assumed a uniform cosmic-ray and gas number density within the target clouds (we 
therefore neglect the temporal, spatial effects within the molecular cloud itself; the whole molec- 
ular clouds becomes instantly a cosmic-ray target). This is a simplification of the model, enough 
however for the aims herein pursued, which is trying to determine the diffusion environment in the 
environment, i.e., between the shell and the cloud, outside the latter. The CR spectrum is given by 

J P (E, r, t) = [cp/4n]f, (3) 

where f(E, r, t) is the distribution function of protons at an instant t and distance r from the 
source. The distribution function satisfies the radial-temporal-energy dependent diffusion equation 
(Ginzburg & Syrovatskii 1964): 

(df/dt) = (D(E)/r 2 )(d/dr)r 2 (df/dr) + (d/dE) (Pf) + Q, (4) 

where P = -dE/dt is the energy loss rate of the particles, Q = Q(E, r, t) is the source function, 
and D(E) is the diffusion coefficient, for which we assume here that it depends only on the par- 
ticle's energy. The energy loss rate are due to ionization and nuclear interactions, with the latter 
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Figure 1. Earlier MAGIC and EGRET (stars and diamonds, respectively), and recent Fermi and VERITAS (squares and upper trianges, respectively) 
measurements of the neighborhood of IC 443 as compared with model predictions for an impulsive and a continuous accelerator, as considered 
in Paper I. The nominal values of parameters for these models are the following: At the MAGIC energy range, the left panel curves show the 
predictions arising from the pp interactions in a cloud of 8000 M G located at 20 (1), 25 (2), and 30 (3) pc, whereas they correspond to 15 (1), 20 
(2), 25 (3), and 30 (4) pc in the right panel. At the EGRET energy range, the curve shows the prediction for a few hundred M Q located at 3^4 pc. 
The right panels show the same results than those in the left, but summing up the different contributions. 



dominating over the former for energies larger than 1 GeV. The nuclear loss rate is P auc = E/t pp , 
with T pp = (n p c k o-p p y l being the timescale for the corresponding nuclear loss, k ~ 0.45 being the 
inelasticity of the interaction, and cr pp being the cross section (Gaisser 1990). Aharonian & Atoyan 
(1996) presented a solution for the diffusion equation for an arbitrary energy loss term, diffusion 
coefficient, and impulsive injection spectrum fi n j(E), such that Q(E, r, t) = N fi^(E)6r6(t). For the 
particular case in which D(E) oc E s and / m j oc E~ a , above ~ 10 GeV, where the cross-section to pp 
interactions is a weak function of E, the general solution is 

f(E, r, t) ~ (N E- a /n 3/2 R 3 dli ) exp [-(a - l)t/r pp - (R/R dit ) 2 } , (5) 

where R d n = 2(D(E)t[exp(tS/T pp ) - l]/[t6/r pp ]) l/2 stands for the radius of the sphere up to which 
the particles of energy E have time to propagate after their injection. In case of continuous injec- 
tion of accelerated particles, given by Q{E, t) = Q E~ a T(t), the previous solution needs to be 
convolved with the function T(t - f) in the time interval < f < t. If the source is described by a 
Heavside function, T(t) = 0(0 Atoyan et al. (1995) have found a general solution for the diffusion 
equation with arbitrary injection spectrum, which with the listed assumptions and for times t less 
than the energy loss time, leads to: 

e' x2 dx. (6) 

r/Rm 

We will assume that a = 2.2 and make use of these solutions in what follows. In Paper I, we gave a 
detailed description of the multi-frequency knowledge on IC 443, impacting on the determination 
of the main parameters entering into the model (e.g., the SNR's age, and molecular environment). 
We refer the reader to that discussion for details. 
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3.1 The model and new data 

3.2 Comparison with nominal models in Paper I 

We start by directly comparing the predictions made in Paper I with the most recent results ob- 
tained by VERITAS and Fermi. In the case of VERITAS, given that their measured SED is com- 
patible with the earlier one obtained by MAGIC, we will no see no significant difference in the 
response of the models. In the case of Fermi, the situation is different because Fermi results ex- 
tended the energy domain of the SED much beyond what was possible for EGRET, and for that 
region the models explored in Paper I, thus, were unconstrained. 

In the models of Paper I, IC 443 was considered both as a continuous accelerator with a rela- 
tivistic proton power of L p = 5 x 10 37 erg s _1 (the proton luminosity is such that the energy injected 
into relativistic CRs through the SNR age is 5 x 10 49 erg), and an impulsive injector with the same 
total power (injection of high energy particles occur in a much shorter time than the SNR age). 
Cosmic-rays were assumed to propagate with a diffusion coefficient at 10 GeV, e.g., D w = 10 26 
cm 2 s _1 , and 6 = 0.5 in a medium of typical density, for which the timescale for nuclear loss t pp is 
orders of magnitude larger than the age of the accelerator. The nominal models explored also took 
assumptions regarding the location (different distances between the SNR shock and the interacting 
clouds were assumed) and molecular mass affected by the cosmic-rays. These assumptions were 
based on the observations of molecular lines towards IC 443 made by, e.g., Cornett et al. (1977), 
De Noyer (1981), Dickman et al. (1992), Seta et al. (1998), and Torres et al. (2003) which conform 
the overall picture: a total mass of ~ 1.1 x 10 4 M Q mostly at the foreground of the remnant, since it 
is found to be absorbing optical and X-ray radiation, with smaller cloud(s) totalizing the remain- 
ing mass located closer to the SNR. Important details are however uncertain, for instance, whether 
there is one or several foreground clouds, the distance between the foreground cloud(s) and the 
SNR shell, the number and specific location(s) of the foreground cloud(s), and their mass distri- 
bution if more than one cloud is there. It was the hoped that Fermi data would elucidate some of 
these parameters, regarding not only the molecular environment but also the diffusion properties 
of the medium, by a posteriori comparison with data. 

Figured] shows the result of the nominal model predictions (theoretical curves are exactly as in 
Paper I, except for the fact we computed and added a contribution of bremsstrahlung radiation (we 
consider this contribution for primary particles with a proton to electron ratio of 150, following 
the standard formulae quoted for instance in Torres 2004) which is only visible at the smallest 
energies in the plot, compared with the newest data. We note that electron bremsstrahlung can 
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hardly explain the whole of the observed IC 443 gamma-ray emission, a conclusion also reached 
by Abdo et al. (2009), and with EGRET data, also by Butt et al. (2003). Since the cross section of 
bremsstrahlung and pion production are similar at the Fermi range, their emission ratio is similar 
to the electron-to-proton ratio. The curves in Figure Q] are based on assuming 8000 M Q at the 
different distances marked in the plot and a few hundred M Q located closer to the SNR (as an 
example, ~700 M Q for the case of an impulsive, and ~300 M Q for a continuous case, located at 
3-4 pc). One can immediately see that what earlier was, particularly in the case of the impulsive 
accelerator, a good agreement between theory and the observations performed by EGRET and 
MAGIC (and also VERITAS) only, is now in disagreement with Fermi data. The spectrum is 
harder than what was suggested by EGRET, presenting an almost flat SED up to 10 GeV, with a 
roll-over in the spectrum between 10 and 100 GeV. Models in Paper I are unable to reproduce the 
details of these trends: In fact, the case of continuous acceleration was already not favored in 
Paper I due to both, the middle age of the remnant and the behavior at the highest energies, which 
were producing a SED much harder than observed and it is now ruled out. We will not consider 
this case any further. In the case of impulsive acceleration, it is at the earlier unexplored region of 
energies, between 10 and 100 GeV, where we find significant deviations between theory and data, 
and there is no model among the ones explored above which can accommodate at the same time 
a SED that is both, sufficiently steep at VHEs to concur with MAGIC/VERITAS observations and 
sufficiently flat one decade earlier in energy to concur with Fermi data. 

3.3 Using Fermi data to constrain model parameters 

What at first sight could seem as a difficult-to- solve failure of the scenario, we find that it is actually 
only the failure of some numerical values of parameters. In particular, differences in the location 
and masses of the overtaken clouds can move the peaks of their corresponding contributions (see 
Aharonian & Atoyan 1996, Gabici et al. 2007, Rodriguez-Marrero et al. 2009 for detailed analysis 
of the dependences). Certainly, kinematic distance estimations are not accurate enough to obtain 
exact separation of the cloud(s) from the SNR shell. Thus, Fermi observations are holding the 
key to make some precisions on the assumptions made in this sense, given that the unknowns can 
affect the final results on the predicted spectra. Using Fermi results we find that a closer (e.g., at 
10 pc) less massive giant cloud (~5300 M©) being overtaken by cosmic-rays diffusing away from 
IC 443 and an smaller amount of molecular material in cloud(s) closer to the SNR shell (e.g., at 4 
pc, with 350 M Q ) produce an excellent match to the whole range of observations, see Figured A 
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Figure 2. As in Figure 1 . Summed results (right) are produced by two main components (shown in the left panel) coming from a giant cloud in 
front of the SNR, which is at least partially overtaken by the diffusing cosmic-rays (~5300 M© at 10 pc) and a closer-to-the-shell cloud, similar 
to the previous examples (in this case, at 4 pc, with 350 M Q ). The dotted (dashed) line at the VHE range corresponds to different normalizations 
(equivalently, to interacting masses of ~4()00 and ~3200 M G at the same distance). The diffusion coefficient is as before, D\q = 10 26 cm 2 s , 

smaller -than the total quoted: 1.1 x 10 4 M Q - amount of mass in the foreground giant molecular 
cloud(s) being overtaken by diffusing cosmic-rays from IC 443 is perfectly possible, given the 
various uncertainties in the absolute position of the cloud, its real number, and the velocity model 
used; and essentially, due to the fact that the total amount of mass correspond to a larger projected 
sky area. Whereas this implies no substantial change to the model, it allows a match with data at 
all high-energy frequencies, see Figure El As in Figure 1, this Figure shows the results produced 
by two main components, one coming from a giant cloud in front of the SNR, which is at least 
partially overtaken by the diffusing cosmic-rays (e.g., ~5300 M at 10 pc) and a closer-to-the-shell 
cloud, of overall magnitudes similar to the previous examples (in this case, at 4 pc, with 350 M Q ). 
The dotted (dashed) line at the VHE range corresponds to different normalizations (equivalently, 
to interacting masses of ~4000 and ~3200 M Q at the same distance). The diffusion coefficient is 
as before, £>i = 10 26 cm 2 s" 1 . 

3.4 Cosmic-ray distributions and their effects 

In Figure [3] we show the distribution of cosmic-rays generated by the impulsive IC 443 at the two 
different distances considered for the molecular mass distribution in one matching model (solid 
black line) of Figure [21 We also plot their ratio with respect to the Earth cosmic-ray distribution. 
It can be seen that the cosmic-ray energy density is greatly enhanced -along the energy range of 
interest — as compared with that in our vicinity, described with an spectrum of the form J e (E) ~ 
2.2£' ( ; 2 y 5 cnT 2 GeV 1 s _1 sr _1 (e.g. Dermer 1986). It can also be seen that significant deviations 
of the cosmic-ray density are obtained in case the diffusion is slower (i.e., D l0 is larger). At a fixed 
SNR age of 30 kyrs, increasing D 10 produces the y-ray emission prediction to displace to smaller 
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Figure 3. Cosmic-ray spectrum generated by the impulsive IC 443 at the two different cloud distances considered in one matching model of Figure 
[2] 10 (solid) and 4 pc (dashed), at the age of the SNR, as a function of energy. Different colors show results for different diffusion coefficient (black, 
D[q = 10 26 cm 2 s" 1 ; and red, D10 = 10 27 cm 2 s" 1 ). The right panel shows the ratio between the cosmic-ray spectra of the top panel, and the 
cosmic-ray spectrum near Earth, as a function of energy. 
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Figure 4. Example of model output with diffusion coefficient scale equal to 10 27 cm 2 s . The different curves represent results for the location of 
giant molecular cloud at 10, 15, 20, 25, 30 pc from the SNR shell, whereas the close-to-the-SNR cloud is at 4 pc. Neither in this nor in any other 
of the models studied varying the parameters with such diffusion coefficient scale, the VHE source spectrum can be reproduced, nor the resulting 
SED in Fermi range is hard enough to match the data. 

energies, typically, until D 10 > D transition , where peaks generated by clouds at large separation (e.g, 
100 pc) displace up and peaks generated by clouds at smaller separation (e.g., 10 pc) displace down 
in the SED (e.g., Rodriguez-Marrero et al. 2009 and references therein). This fact implies that for 
the range of distances to the giant and close-to-the-SNR molecular clouds being considered (10-30 
pc, and 2-6 pc respectively) there is no solution with large D l0 able to fit the whole range of data. 
This was already hinted at in Paper I, where just using MAGIC data we found that it was possible 
to put an strong constraint over the diffusion timescale: D 10 should be of the order of 10 26 cm 2 s _1 
since if the separation between the giant cloud and the SNR is >10 pc, an slower diffusion would 
not allow sufficient high energy particles to reach the target material and it would be impossible 
to reproduce the VHE data. On the other hand, given that there is a displacement between the 
centroid positions of EGRET/Fermi and VHE sources and that molecular material is absorbing 
lower frequency emission from the remnant, the separation between the foreground cloud(s) and 
SNR shell can not be much smaller than 10 pc. The current Fermi data emphasizes this conclusion. 
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Figure 5. Contour plot depicting the position of the peak of the SED generated by a 30000 years old injection interacting with clouds at different 
distances, for a range of diffusion coefficient scale, Dio 

Figure |4] shows an example of a full range of models we constructed with D\q = 10 27 cm 2 s -1 , and 
its disagreement with data. Note that even if rescaling the curves assuming e.g., a much higher 
molecular mass (which would in itself be in conflict with multi-frequency observations), it is not 
possible to obtain a good fit across the whole range of observations. 

Figure [5] shows, as contour plots, the energy at which the maximum of the SED is found for 
the cases of impulsive acceleration of cosmic -rays. The age corresponds to the SNR like IC 443, 
and the cosmic-rays are interacting with clouds at different distances, for a range of diffusion 
coefficient scale, D w . This plot is useful to notice which is the solution we are finding and its 
degree of uncertainty/degeneracy: in order to fit the combined MAGIC/VERITAS and Fermi data 
we would need a giant cloud producing a peak at about the Fermi spectral turnover (what means, 
looking at the plot, either a very large separation between the cloud and the SNR shell for a high 
Dio, (what we discarded because it is not possible to fit the VHE MAGIC/VERITAS data in this 
configuration), or a smaller distance with a faster diffusion, i.e. a lower D 10 , which is the solution 
we can still promote. 

3.5 More on degeneracies and uncertainties in parameter estimation 

Figure [6] explores the range of parameters around the solutions matching the observational data; 
giving a feeling of the degeneracies (or uncertainties) within which this model provides a reason- 
able agreement with observations. The values of masses and diffusion coefficients used in Figure 
[6] to obtain good data-matching given the distances to each of the clouds are given in Table 1. Fits 
could be considered good for D GM c between 9 and 11 pc. For an average distance of 10 pc, the 
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Figure 6. Examples of solutions around the main values discussed, exploring the degeneracies (or uncertainties) in determining the numerical 
values of model parameters matching the observational data. The order of the panels in this plot, top left to bottom right, corresponds with the 
parameters described in Table 1 . 

mass in the close-to-the-SNR cloud (or clouds) decreases the farthest the latter is. For these cases, 
good solutions with D GMC between 9 and 1 1 pc can always be found adjusting other parameters. 
Our average model explored in Figure [2] corresponds to D GMC = 10 pc, d snr =A pc, with the three 
curves constructed with -5300, -4000, -3200 M©, and M snr = 350 M G . The results in Figure [6] 
and Table 1 show that the smaller the diffusion coefficient, the fit at VHEs worsens, overpredicting 
the data. Correcting this via a mass adjustment, would in turn make for a poor fit at lower ener- 
gies; what in practice imposes a lower limit to D w . On the other hand, when D l0 increases the 
VHE spectra is quickly underpredicted, and again, correcting this via a mass adjustment would 
in turn make for a poor fit at lower energies. In summary, in order for this model to match the 
multi-frequency observational data, the range of variation in the parameters gets constrained as 
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Table 1. Main model parameters for solutions shown in Figure[6] Dqmc an d d mr are the distance to the GMC and the closer-to-the-SNR molecular 
clouds. The three / values quoted define Mqmc — (1//) 8000 M Q . The three groups explore different degeneracies: in position of the GMC, of the 
smaller cloud, and on the diffusion coefficient. Model 3 is not shown in Figure|6]but rather in Figure[2] 
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9 < D GM c ^11 P c ; 3 < d snr < 6 pc, and D w ~ 10 26 cm 2 s _1 ; thus constituting a direct estimation 
of D l0 and the molecular environment in the IC 443 vicinity under the assumed validity of this 
model. We have mentioned above that we assumed the y-ray emissivity was constant within the 
clouds; i.e. we are assuming that there is no significant cosmic-ray gradient in the target. This as- 
sumption is an approximation, which is better when the size of the cloud is less than the distance 
to the accelerator and the diffusion coefficients inside and outside the cloud are not significantly 
different (or even if they are, the proton-proton timescale is larger than the time it takes for cosmic 
rays to overtake the whole cloud). In the case of IC 443, these conditions can be accommodated 
for the solutions in Table 1, except perhaps for the very massive cloud located close to the SNR at 
d snr = 2 pc; for which would imply a cloud average density higher than usually found, although 
even this might also be possible given the small scale clumps found therein (e.g., see Rosado et al. 
2007). 

3.6 Influence of the 5-parameter 

We have also made an exploration of other parameters of the model, as for instance, those influenc- 
ing the way in which the diffusion coefficient varies with energy (the parameter 5), or the injection 
spectrum of cosmic rays (referred to as a), and came to the conclusion that their corresponding 
values are rather constrained. For instance, the 8 parameter is expected to be around 6 = 0.4 - 0.7 
(e.g., Berezinskii et a. 1990) and a typical value of 0.5 is usually assumed. Figure [7] gives account 
of how small variations in 5 change the slope of good-fitting solutions to the high and very-high 
energy data. One can see that for steeper ^-parameters, of course, steeper y-ray spectrum are found. 
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Figure 7. Comparing y-ray yields with different S parameters, from left to right 6 = 0.4, 0.5, 0.6, and 0.7. Other parameters are as in Figure 2. 

If the masses of the molecular clouds are maintained and 6 is larger, in order to have a good fit one 
would need an even lower D w , lower than D 10 = 10 26 cm 2 s" 1 (see Figure©, making the solution 
less feasible. 



3.7 Uncertainties due to the cross section parameterization 

We have checked whether changes in the cross section parameterization can produce significant 
variance in the results. In the appendix of Domingo Santamana and Torres (2005), the different 
predicted yields in y-rays obtained when using alternate cross section parameterizations known by 
then were compared among themselves and with data. The parameterizations therein considered 
were Kamae et al.'s (2005); the ^-functional form by Aharonian & Atoyan (1996) that is used 
above, Stephen and Badwhar's (1981), and Blattnig et al.'s (2000a,b). It was found that Kamae's 
and the ^-functional form were very close to each other, as seen in Fig. 1 1 of that paper, which 
showed the y-ray emissivities obtained with the corresponding use of each of the parameterizations 
of the cross section. In that paper, it was also found that neither Stephen and Badwhar's (1981) 
nor Blattnig et al.'s (2000a,b) were appropriate for their use in broad-band high-energy modeling 
such as the one we pursue here. More recently, Kelner et al. (2006) presented a new approach for 
obtaining the cross section in pp interactions. These authors used 2 shapes for representing the 
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Figure 8. Comparison (left) and ratio (right) of the cross section parameterizations used in the figures above, the ^-functional form by Aharonian 
& Atoyan (1996), with that of Kelner et al. (2006). 



cross section, separated in energy. At low energies (up to 100 GeV), Kelner et al. approach uses a 
slightly modified but similarly-shaped ^-functional form. At high-energy, its approach is different, 
and consists of presenting an analytical shape fitting of the results of the simulations of energy 
distribution of n mesons by the SYBILL code. Figure [8] shows a comparison of the cross section 
parameterizations used in the figures above, the ^-functional form by Aharonian & Atoyan (1996), 
with that of Kelner et al. (2006). Differences are within 20%, with the 5-functional approximation 
being larger. 

The concomitant change in flux predictions, produced only because of a different use of a cross 
section parameterization, then, can be reabsorbed as part of the uncertainty in the determination 
of the model. For instance, Figure [9] shows its impact in two ways: The left panel shows several 
alternatives for the position of the large (TeV-producing) cloud in the model; located at 10, 15, and 
20 pc (here we maintain the mass of this cloud fixed and look only at the shape of the curves for 
different distances). It is clear that the change in cross section parameterization does not make any 
of the previously unfeasible models, feasible, and again single out a distance of about 10 pc from 
the SNR shell to the giant TeV-producing cloud for obtaining a good fit. The right panel assumes 
this distance of 10 pc and explores the uncertainty in the determination of the cloud mass. The 
parameters therein shown are 4 pc, and 350 M for the close-to-the-remnant cloud, i.e., the same 
as above, and 10 pc and 7272, 5333, 4210 M Q for the TeV-producing giant cloud. 



3.8 Computation of secondaries other than photons 

With the use of the Kelner et al. (2006) parameterization one can also readily compute secondaries 
other than photons, and this is shown in Figure \10\ Gabici et al. (2009) showed that secondary 
electrons produced within clouds of a wide range of parameters can escape without being affected 
by significant losses; i.e. that the propagation time through the cloud for cosmic-ray electrons is 
shorter than the energy loss time for particles energies between ~ 100 MeV and few hundreds 
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Figure 9. Left: y-ray flux results using the Kelner et al. (2006) approximation for different distances from the shell to the TeV-producing cloud, 10 
(solid), 15 (dotted) and 20 (dashed) pc. The close-to-the-remnant cloud is fixed at 4 pc and contains 350 M G - changes in these latter value do not 
improve the overall fit. Right: For the best-fitting distance models, y-ray flux results using the Kelner et al. (2006) for different values of the giant 
cloud mass (see text for details). 
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Figure 10. Electrons (electrons and positrons are shown together), photons and two flavors of neutrinos produced within the clouds considered 
nearby IC 443, using a set of parameters shown in Figure|9] right panel, with mass of the giant cloud equal to 7272 M Q . The and v e neutrino 
curves show the particle and the anti-particle flux together. Data should only be compared with the photon curve. 



TeV. There would be, then, little effect of the secondary electrons produced on the non-thermal 
emission from the cloud. In addition, for typical densities of clouds, in the several hundreds to 
several thousands particles per cm 3 , the dominant energy loss from ~ 100 MeV and ~ 10 TeV, 
would be bremsstrahlung and not synchrotron. 

A conclusive proof of the hadronic nature of the gamma-ray emission could however come 
from the detection of neutrinos. Neutrino telescopes search for up-going muons produced deep in 
the Earth, and are mainly sensitive to the incoming flux of v M and v M . The finished ICECUBE, for 
example, will consist of 4800 photomultipliers, arranged on 80 strings placed at depths between 
1400 and 2400 m under the South Pole ice (e.g., Halzen 2006). The strings will be located in 
a regular space grid covering a surface area of 1 km 2 . Each string will have 60 optical modules 
(OM) spaced 17 m apart. The number of OMs which have seen at least one photon (from Cerenkov 
radiation produced by the muon which resulted from the interaction of the incoming v in the earth 
and ice crust) is called the channel multiplicity, N ch . The multiplicity threshold is set to N c h = 10, 
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which corresponds to an energy threshold of 200 GeV. The angular resolution of ICECUBE will 
be around ~ 0.7°. 

A first estimation of the event rate of the atmospheric v-background that will be detected in the 
search bin can be obtained as (e.g., Anchordoqui et al. 2003) 



where A eff is the effective area of the detector, AQ. w 1.5 x 10^ 4 sr is the angular size of the 
search bin, and d^/dEy < 0.2 (E v /GeW)~ %21 GeV -1 cnr 2 s _1 sr -1 is the + v M atmospheric 
v-flux (Volkova 1980, Lipari 1993). Here, P V ^(E V ) denotes the probability that a v of energy E v 
on a trajectory through the detector, produces a muon. For E v ~ 1 — 10 3 GeV, this probability is 
* 3.3 x 10~ 13 (£ v /GeV) 2 - 2 , whereas for E v > 1 TeV, P V ^{E V ) » 1.3 X 10~ 6 (£ y /TeV) 8 (Gaisser 
et al. 1995). On the other hand, the v-signal is similarly obtained as 



where (F v + F v ) is the incoming v^-flux. In the previous integrals we use both expressions for 
Py^(E v ) according to the energy, and integrate from 200 GeV up to 10 TeV. The effect of v- 
oscillations is taken into account following table 2 of Cavassini et al. (2006), where the oscillation 
probability in the average vacuum oscillation hypothesis is given. It is is assumed that the inter- 
conversion probability between flavors and between anti-flavors is the same. As an effect of oscil- 
lations, the flavor composition of all the expected fluxes for each flavor are within 50% of each 
other. Using the former formulae and the secondary computation shown in Figure [TO] we find that 
the number of muon neutrino signal events is 0.6 per year of observation, still significantly below 
than the estimation of the number of background events, which under the previous provisions is 
6.4 along the same period, with the full ICECUBE array. If we consider only events above 1 TeV, 
the expected signal is 0.25 year 1 , and the computed background is 1.92 year 1 . ICECUBE does 
not seem to be able to distinguish this signal in reasonable integration times, at least within the 
reach of this simplified treatment of the detector. 

4 CONCLUDING REMARKS 

The recent observations of the IC 443 environment made by AGILE, Fermi, and VERITAS at the 
GeV and TeV energies are spectrally consistent with the interpretation of cosmic-ray interactions 
with a giant molecular cloud lying in front of the remnant. This scenario would be producing 
no significant counterpart at lower energies at that spot, and would then be leading to a natural 




(7) 




(8) 
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interpretation of the dislocation between the centroids of the detections at the different energy 
bands. Use of the latest data allowed to estimate, within the assumed validity and framework of 
this model, the diffusion characteristics in this environment, showing that the diffusion coefficient 
is lower; the cosmic-ray density is higher, than the Earth- values of these magnitudes. Uncertain- 
ties in amount and localization of target molecular mass still remains as does also in the density 
at which this molecular material is found (e.g., the uncertainty in the cosmic-ray-overtaken mass 
discussed above is about 100% for matching models at the extremes of this parameter). But even 
allowing for a range this large, the model could accommodate some but not all variations in other 
parameters, with the values of D w and distances from the SNR shell seemingly being solid con- 
straints. 
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